Подготовка к работе

Нам понадобятся следующие пакеты.

2. Подготовим данные

Проверим отсутсвующие значения и выбросим их при необходимости.

numberOfNA <- length(which(is.na(Boston)==T))
if(numberOfNA>0) {
  Boston <- Boston[complete.cases(Boston),]
}

Создадим данные для тестирования модели.

set.seed(123)
split <- sample.split(Boston,SplitRatio = 0.75) #assigns booleans to a new coloumn based on the split ratio
train <- subset(Boston,split==TRUE)
test <- subset(Boston,split==FALSE)

3. Разведочный анализ.

Посмотрим структуру датасета и базовую статистику.

str(Boston)

## 'data.frame':    506 obs. of  14 variables:
##  $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
##  $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
##  $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
##  $ chas   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
##  $ rm     : num  6.58 6.42 7.18 7 7.15 ...
##  $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
##  $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
##  $ rad    : int  1 2 2 3 3 3 5 5 5 5 ...
##  $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
##  $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
##  $ black  : num  397 397 393 395 397 ...
##  $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
##  $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...

summary(Boston)

##       crim                zn             indus            chas        
##  Min.   : 0.00632   Min.   :  0.00   Min.   : 0.46   Min.   :0.00000  
##  1st Qu.: 0.08204   1st Qu.:  0.00   1st Qu.: 5.19   1st Qu.:0.00000  
##  Median : 0.25651   Median :  0.00   Median : 9.69   Median :0.00000  
##  Mean   : 3.61352   Mean   : 11.36   Mean   :11.14   Mean   :0.06917  
##  3rd Qu.: 3.67708   3rd Qu.: 12.50   3rd Qu.:18.10   3rd Qu.:0.00000  
##  Max.   :88.97620   Max.   :100.00   Max.   :27.74   Max.   :1.00000  
##       nox               rm             age              dis        
##  Min.   :0.3850   Min.   :3.561   Min.   :  2.90   Min.   : 1.130  
##  1st Qu.:0.4490   1st Qu.:5.886   1st Qu.: 45.02   1st Qu.: 2.100  
##  Median :0.5380   Median :6.208   Median : 77.50   Median : 3.207  
##  Mean   :0.5547   Mean   :6.285   Mean   : 68.57   Mean   : 3.795  
##  3rd Qu.:0.6240   3rd Qu.:6.623   3rd Qu.: 94.08   3rd Qu.: 5.188  
##  Max.   :0.8710   Max.   :8.780   Max.   :100.00   Max.   :12.127  
##       rad              tax           ptratio          black       
##  Min.   : 1.000   Min.   :187.0   Min.   :12.60   Min.   :  0.32  
##  1st Qu.: 4.000   1st Qu.:279.0   1st Qu.:17.40   1st Qu.:375.38  
##  Median : 5.000   Median :330.0   Median :19.05   Median :391.44  
##  Mean   : 9.549   Mean   :408.2   Mean   :18.46   Mean   :356.67  
##  3rd Qu.:24.000   3rd Qu.:666.0   3rd Qu.:20.20   3rd Qu.:396.23  
##  Max.   :24.000   Max.   :711.0   Max.   :22.00   Max.   :396.90  
##      lstat            medv      
##  Min.   : 1.73   Min.   : 5.00  
##  1st Qu.: 6.95   1st Qu.:17.02  
##  Median :11.36   Median :21.20  
##  Mean   :12.65   Mean   :22.53  
##  3rd Qu.:16.95   3rd Qu.:25.00  
##  Max.   :37.97   Max.   :50.00

Здесь мы можем видеть, что переменные «crim» и «black» принимают широкий диапазон значений.

Переменные «crim», «zn», «rm» и «black» имеют большую разницу между их медианной и средним значением, что указывает на множество выбросов в соответствующих переменных.

par(mfrow = c(1, 4))
boxplot(Boston$crim, main='crim',col='Sky Blue')
boxplot(Boston$zn, main='zn',col='Sky Blue')
boxplot(Boston$rm, main='rm',col='Sky Blue')
boxplot(Boston$black, main='black',col='Sky Blue')

Как и было предсказано, мы видим много выбросов.

Посчитаем корреляцию между признаками.

corr_matrix <- cor(Boston)
corrplot.mixed(corr_matrix)

Из корреляционной матрицы можно сделать следующие наблюдения:

Средняя стоимость домов, занимаемых владельцами, увеличивается по мере увеличения среднего количества комнат на жилище и уменьшается, если процент населения с более низким статусом в этом районе увеличивается.
Концентрация NOx или оксидов азота (ppm) возрастает с увеличением доли акров, не относящихся к розничной торговле, на город и доли жилых единиц, построенных до 1940 года.
rad и tax имеют сильную положительную корреляцию 0,91, что означает, что по мере увеличения доступности радиальных автомагистралей, полная ставка налога на имущество из расчета на $ 10 000 также увеличивается.
crim тесно связан с переменными rad и tax, что означает, что по мере увеличения доступности радиальных автомагистралей увеличивается уровень преступности на душу населения.
indus имеет сильную положительную корреляцию с NOx, что подтверждает мнение о высокой концентрации оксидов азота в промышленных зонах.

Визуализируем данные.

Boston %>%
  gather(key, val, -medv) %>%
  ggplot(aes(x = val, y = medv)) +
  geom_point() +
  stat_smooth(method = "lm", se = TRUE, col = "blue") +
  facet_wrap(~key, scales = "free") +
  theme_gray() +
  ggtitle("Scatter plot of dependent variables vs Median Value (medv)")

## `geom_smooth()` using formula 'y ~ x'

4. Построение модели, прогнозирующей стоимость дома

Посмотрим распределение переменной medv.

qplot(x=medv,data=Boston,geom='histogram')

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

На данной гистограмме распределения средней цены дома мы видим, что данные немного смещены влево. Однако, линейную регрессию нельзя строить на ненормально распределенных данных. Чтобы исправить распределние прологарифмируем параметр medv.

qplot(x=log(medv),data=Boston,geom='histogram')

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

После логарифмирования у нас остался небольшой сдвиг вправо, но распределение выглядит более “нормальным”.

Задание 1. Построение полной модели.

Стандартизируем предикторы:

# сделаем наличие реки фактором 
Boston$chas <- as.factor(Boston$chas)
stand_Boston <- Boston %>% mutate_at(c('crim', 'zn', 'indus', 'nox', 'age', 'dis', 'tax', 'ptratio', 'black', 'lstat'), scale)

model1_stand = lm(log(medv)~., data= stand_Boston)
summary(model1_stand)

## 
## Call:
## lm(formula = log(medv) ~ ., data = stand_Boston)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.73361 -0.09747 -0.01657  0.09629  0.86435 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.320438   0.104407  22.225  < 2e-16 ***
## crim        -0.088351   0.011315  -7.808 3.52e-14 ***
## zn           0.027345   0.012815   2.134 0.033349 *  
## indus        0.016923   0.016886   1.002 0.316755    
## chas1        0.100888   0.034486   2.925 0.003598 ** 
## nox         -0.090199   0.017717  -5.091 5.07e-07 ***
## rm           0.090833   0.016728   5.430 8.87e-08 ***
## age          0.005928   0.014883   0.398 0.690567    
## dis         -0.103364   0.016811  -6.149 1.62e-09 ***
## rad          0.014267   0.002656   5.373 1.20e-07 ***
## tax         -0.105466   0.025368  -4.157 3.80e-05 ***
## ptratio     -0.082856   0.011337  -7.309 1.10e-12 ***
## black        0.037757   0.009815   3.847 0.000135 ***
## lstat       -0.207344   0.014496 -14.304  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1899 on 492 degrees of freedom
## Multiple R-squared:  0.7896, Adjusted R-squared:  0.7841 
## F-statistic: 142.1 on 13 and 492 DF,  p-value: < 2.2e-16

par(mfrow=c(2,2))
plot(model1_stand)

Полная модель:

R-squared составляет 0.7841
F-statistic составляет 142.1

Полная модель объясняет 78% изменчивости. Протестируем, есть ли параметры, которые незначительно влияют на модель. Для этого воспользуемся критерием Акаике.

step = stepAIC(model1_stand,direction = 'both')

## Start:  AIC=-1667.19
## log(medv) ~ crim + zn + indus + chas + nox + rm + age + dis + 
##     rad + tax + ptratio + black + lstat
## 
##           Df Sum of Sq    RSS     AIC
## - age      1    0.0057 17.755 -1669.0
## - indus    1    0.0362 17.786 -1668.2
## <none>                 17.749 -1667.2
## - zn       1    0.1643 17.914 -1664.5
## - chas     1    0.3088 18.058 -1660.5
## - black    1    0.5339 18.283 -1654.2
## - tax      1    0.6235 18.373 -1651.7
## - nox      1    0.9351 18.684 -1643.2
## - rad      1    1.0413 18.791 -1640.3
## - rm       1    1.0637 18.813 -1639.7
## - dis      1    1.3639 19.113 -1631.7
## - ptratio  1    1.9270 19.676 -1617.0
## - crim     1    2.1995 19.949 -1610.1
## - lstat    1    7.3809 25.130 -1493.2
## 
## Step:  AIC=-1669.03
## log(medv) ~ crim + zn + indus + chas + nox + rm + dis + rad + 
##     tax + ptratio + black + lstat
## 
##           Df Sum of Sq    RSS     AIC
## - indus    1    0.0363 17.791 -1670.0
## <none>                 17.755 -1669.0
## + age      1    0.0057 17.749 -1667.2
## - zn       1    0.1593 17.914 -1666.5
## - chas     1    0.3138 18.069 -1662.2
## - black    1    0.5431 18.298 -1655.8
## - tax      1    0.6205 18.376 -1653.7
## - nox      1    0.9645 18.720 -1644.3
## - rad      1    1.0356 18.791 -1642.3
## - rm       1    1.1452 18.900 -1639.4
## - dis      1    1.5471 19.302 -1628.8
## - ptratio  1    1.9224 19.677 -1619.0
## - crim     1    2.1988 19.954 -1612.0
## - lstat    1    8.1949 25.950 -1479.0
## 
## Step:  AIC=-1670
## log(medv) ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + 
##     black + lstat
## 
##           Df Sum of Sq    RSS     AIC
## <none>                 17.791 -1670.0
## + indus    1    0.0363 17.755 -1669.0
## + age      1    0.0058 17.786 -1668.2
## - zn       1    0.1451 17.936 -1667.9
## - chas     1    0.3399 18.131 -1662.4
## - black    1    0.5344 18.326 -1657.0
## - tax      1    0.6139 18.405 -1654.8
## - nox      1    0.9350 18.726 -1646.1
## - rad      1    1.0088 18.800 -1644.1
## - rm       1    1.1171 18.909 -1641.2
## - dis      1    1.7385 19.530 -1624.8
## - ptratio  1    1.8862 19.678 -1621.0
## - crim     1    2.2229 20.014 -1612.4
## - lstat    1    8.1604 25.952 -1481.0

step$anova

## Stepwise Model Path 
## Analysis of Deviance Table
## 
## Initial Model:
## log(medv) ~ crim + zn + indus + chas + nox + rm + age + dis + 
##     rad + tax + ptratio + black + lstat
## 
## Final Model:
## log(medv) ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + 
##     black + lstat
## 
## 
##      Step Df    Deviance Resid. Df Resid. Dev       AIC
## 1                              492   17.74938 -1667.194
## 2   - age  1 0.005723781       493   17.75510 -1669.031
## 3 - indus  1 0.036264380       494   17.79137 -1669.999

step_model = lm(log(medv) ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + black + lstat,data=stand_Boston)

summary(step_model)

## 
## Call:
## lm(formula = log(medv) ~ crim + zn + chas + nox + rm + dis + 
##     rad + tax + ptratio + black + lstat, data = stand_Boston)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.73400 -0.09460 -0.01771  0.09782  0.86290 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.328996   0.101070  23.043  < 2e-16 ***
## crim        -0.088757   0.011298  -7.856 2.49e-14 ***
## zn           0.025361   0.012637   2.007 0.045308 *  
## chas1        0.105148   0.034228   3.072 0.002244 ** 
## nox         -0.083634   0.016414  -5.095 4.97e-07 ***
## rm           0.090673   0.016281   5.569 4.20e-08 ***
## dis         -0.108878   0.015671  -6.948 1.18e-11 ***
## rad          0.013446   0.002540   5.293 1.82e-07 ***
## tax         -0.094023   0.022774  -4.129 4.28e-05 ***
## ptratio     -0.081025   0.011196  -7.237 1.77e-12 ***
## black        0.037679   0.009781   3.852 0.000133 ***
## lstat       -0.204262   0.013570 -15.053  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1898 on 494 degrees of freedom
## Multiple R-squared:  0.7891, Adjusted R-squared:  0.7844 
## F-statistic: 168.1 on 11 and 494 DF,  p-value: < 2.2e-16

Новая модель:

R-squared составляет 0.7864
F-statistic составляет 104.3

С помощью функции stepAIC мы последовательно рассчитали AIC для некоторых моделей и удалили предикторы, которые незначительно влияют на модель, что позволило несколько увеличить R-squared.

Задание 2. Проведем диагностику полной модели.

Проверим на мультиколлениарность, которую мы заметили ранее.

library(car)
vif(model1_stand)

##     crim       zn    indus     chas      nox       rm      age      dis 
## 1.792192 2.298758 3.991596 1.073995 4.393720 1.933744 3.100826 3.955945 
##      rad      tax  ptratio    black    lstat 
## 7.484496 9.008554 1.799084 1.348521 2.941491

Для дальнейшей работы нужно удалить предикторы с VIF > 5 (rad, tax).

Проверим остатки:

plot(model1_stand$residuals)
abline(h=0,col='red')

Мы видим что остатки распределены ненормально.

Найдем выбросы:

outlierTest(model1_stand)

##      rstudent unadjusted p-value Bonferroni p
## 413  4.781438         2.3041e-06    0.0011659
## 372  4.308070         1.9900e-05    0.0100700
## 373  4.016731         6.8256e-05    0.0345370
## 402 -3.946653         9.0820e-05    0.0459550

Мы видим 4 выброса, влияющие на нашу модель.

Наибольший вклад вносит переменная lstat - -0.207344.

Задание 3. Постройте график предсказаний стоимости от переменной, которая обладает

наибольшим по модулю коэффициентом.

mod_hi <- lm(medv ~ as.numeric(lstat), data = stand_Boston)

test$predicted.medv <- predict(mod_hi,test)
pl1 <-test %>% 
  ggplot(aes(medv,predicted.medv)) +
  geom_point(alpha=0.5) + 
  stat_smooth(aes(colour='black')) +
  xlab('Значения medv') +
  ylab('Предсказанные значения medv') +
  theme_bw()

ggplotly(pl1)

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Дополнительное задание.

Продолжим преобразовывать модель на НЕстандартизованных данных.

# сделаем полную модель на НЕстандартизованных данных
model1 = lm(log(medv)~., data= Boston)
summary(model1_stand)

## 
## Call:
## lm(formula = log(medv) ~ ., data = stand_Boston)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.73361 -0.09747 -0.01657  0.09629  0.86435 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.320438   0.104407  22.225  < 2e-16 ***
## crim        -0.088351   0.011315  -7.808 3.52e-14 ***
## zn           0.027345   0.012815   2.134 0.033349 *  
## indus        0.016923   0.016886   1.002 0.316755    
## chas1        0.100888   0.034486   2.925 0.003598 ** 
## nox         -0.090199   0.017717  -5.091 5.07e-07 ***
## rm           0.090833   0.016728   5.430 8.87e-08 ***
## age          0.005928   0.014883   0.398 0.690567    
## dis         -0.103364   0.016811  -6.149 1.62e-09 ***
## rad          0.014267   0.002656   5.373 1.20e-07 ***
## tax         -0.105466   0.025368  -4.157 3.80e-05 ***
## ptratio     -0.082856   0.011337  -7.309 1.10e-12 ***
## black        0.037757   0.009815   3.847 0.000135 ***
## lstat       -0.207344   0.014496 -14.304  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1899 on 492 degrees of freedom
## Multiple R-squared:  0.7896, Adjusted R-squared:  0.7841 
## F-statistic: 142.1 on 13 and 492 DF,  p-value: < 2.2e-16

par(mfrow=c(2,2))
plot(model1)

Новая модель:

R-squared составляет 0.7841
F-statistic составляет 142.1

step = stepAIC(model1,direction = 'both')

## Start:  AIC=-1667.19
## log(medv) ~ crim + zn + indus + chas + nox + rm + age + dis + 
##     rad + tax + ptratio + black + lstat
## 
##           Df Sum of Sq    RSS     AIC
## - age      1    0.0057 17.755 -1669.0
## - indus    1    0.0362 17.786 -1668.2
## <none>                 17.749 -1667.2
## - zn       1    0.1643 17.914 -1664.5
## - chas     1    0.3088 18.058 -1660.5
## - black    1    0.5339 18.283 -1654.2
## - tax      1    0.6235 18.373 -1651.7
## - nox      1    0.9351 18.684 -1643.2
## - rad      1    1.0413 18.791 -1640.3
## - rm       1    1.0637 18.813 -1639.7
## - dis      1    1.3639 19.113 -1631.7
## - ptratio  1    1.9270 19.676 -1617.0
## - crim     1    2.1995 19.949 -1610.1
## - lstat    1    7.3809 25.130 -1493.2
## 
## Step:  AIC=-1669.03
## log(medv) ~ crim + zn + indus + chas + nox + rm + dis + rad + 
##     tax + ptratio + black + lstat
## 
##           Df Sum of Sq    RSS     AIC
## - indus    1    0.0363 17.791 -1670.0
## <none>                 17.755 -1669.0
## + age      1    0.0057 17.749 -1667.2
## - zn       1    0.1593 17.914 -1666.5
## - chas     1    0.3138 18.069 -1662.2
## - black    1    0.5431 18.298 -1655.8
## - tax      1    0.6205 18.376 -1653.7
## - nox      1    0.9645 18.720 -1644.3
## - rad      1    1.0356 18.791 -1642.3
## - rm       1    1.1452 18.900 -1639.4
## - dis      1    1.5471 19.302 -1628.8
## - ptratio  1    1.9224 19.677 -1619.0
## - crim     1    2.1988 19.954 -1612.0
## - lstat    1    8.1949 25.950 -1479.0
## 
## Step:  AIC=-1670
## log(medv) ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + 
##     black + lstat
## 
##           Df Sum of Sq    RSS     AIC
## <none>                 17.791 -1670.0
## + indus    1    0.0363 17.755 -1669.0
## + age      1    0.0058 17.786 -1668.2
## - zn       1    0.1451 17.936 -1667.9
## - chas     1    0.3399 18.131 -1662.4
## - black    1    0.5344 18.326 -1657.0
## - tax      1    0.6139 18.405 -1654.8
## - nox      1    0.9350 18.726 -1646.1
## - rad      1    1.0088 18.800 -1644.1
## - rm       1    1.1171 18.909 -1641.2
## - dis      1    1.7385 19.530 -1624.8
## - ptratio  1    1.8862 19.678 -1621.0
## - crim     1    2.2229 20.014 -1612.4
## - lstat    1    8.1604 25.952 -1481.0

step$anova

## Stepwise Model Path 
## Analysis of Deviance Table
## 
## Initial Model:
## log(medv) ~ crim + zn + indus + chas + nox + rm + age + dis + 
##     rad + tax + ptratio + black + lstat
## 
## Final Model:
## log(medv) ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + 
##     black + lstat
## 
## 
##      Step Df    Deviance Resid. Df Resid. Dev       AIC
## 1                              492   17.74938 -1667.194
## 2   - age  1 0.005723781       493   17.75510 -1669.031
## 3 - indus  1 0.036264380       494   17.79137 -1669.999

step_model = lm(log(medv) ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + black + lstat,data=Boston)

summary(step_model)

## 
## Call:
## lm(formula = log(medv) ~ crim + zn + chas + nox + rm + dis + 
##     rad + tax + ptratio + black + lstat, data = Boston)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.73400 -0.09460 -0.01771  0.09782  0.86290 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.0836823  0.2030491  20.112  < 2e-16 ***
## crim        -0.0103187  0.0013134  -7.856 2.49e-14 ***
## zn           0.0010874  0.0005418   2.007 0.045308 *  
## chas1        0.1051484  0.0342285   3.072 0.002244 ** 
## nox         -0.7217440  0.1416535  -5.095 4.97e-07 ***
## rm           0.0906728  0.0162807   5.569 4.20e-08 ***
## dis         -0.0517059  0.0074420  -6.948 1.18e-11 ***
## rad          0.0134457  0.0025405   5.293 1.82e-07 ***
## tax         -0.0005579  0.0001351  -4.129 4.28e-05 ***
## ptratio     -0.0374259  0.0051715  -7.237 1.77e-12 ***
## black        0.0004127  0.0001071   3.852 0.000133 ***
## lstat       -0.0286039  0.0019002 -15.053  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1898 on 494 degrees of freedom
## Multiple R-squared:  0.7891, Adjusted R-squared:  0.7844 
## F-statistic: 168.1 on 11 and 494 DF,  p-value: < 2.2e-16

Новая модель:

R-squared составляет 0.7844
F-statistic составляет 168.1

library(car)
vif(step_model)

##     crim       zn     chas      nox       rm      dis      rad      tax 
## 1.789704 2.239229 1.059819 3.778011 1.834806 3.443420 6.861126 7.272386 
##  ptratio    black    lstat 
## 1.757681 1.341559 2.581984

outlierTest(step_model)

##      rstudent unadjusted p-value Bonferroni p
## 413  4.775740         2.3647e-06    0.0011965
## 372  4.340845         1.7232e-05    0.0087195
## 373  4.007183         7.0941e-05    0.0358960
## 402 -3.952037         8.8814e-05    0.0449400

У нас есть 4 выброса, которые мы удалим.

Boston_1 = Boston[-c(413,372, 373, 402),,drop=T]# drop the points
row.names(Boston_1)=1:nrow(Boston_1)

model2 = lm(log(medv)~crim+zn+chas+nox+dis+ptratio+black+lstat,data=Boston)
summary(model2)

## 
## Call:
## lm(formula = log(medv) ~ crim + zn + chas + nox + dis + ptratio + 
##     black + lstat, data = Boston)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.70805 -0.12147 -0.01373  0.10405  0.82894 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.7010244  0.1421502  33.071  < 2e-16 ***
## crim        -0.0080105  0.0012719  -6.298 6.65e-10 ***
## zn           0.0013826  0.0005565   2.484 0.013311 *  
## chas1        0.1302984  0.0363485   3.585 0.000371 ***
## nox         -0.7322562  0.1352610  -5.414 9.63e-08 ***
## dis         -0.0578514  0.0077684  -7.447 4.25e-13 ***
## ptratio     -0.0376285  0.0048687  -7.729 6.06e-14 ***
## black        0.0002958  0.0001120   2.642 0.008505 ** 
## lstat       -0.0353753  0.0017219 -20.545  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2025 on 497 degrees of freedom
## Multiple R-squared:  0.7584, Adjusted R-squared:  0.7545 
## F-statistic:   195 on 8 and 497 DF,  p-value: < 2.2e-16

plot(model2)

Новая модель:

R-squared составляет 0.7545
F-statistic составляет 195

Не смотря на то, что мы устранили мульитколлинеанрость и незначимые предикторы, предсказания модели стали немного хуже (75%).

vif(model2)

##     crim       zn     chas      nox      dis  ptratio    black    lstat 
## 1.473811 2.074443 1.049493 3.024847 3.294717 1.367971 1.286382 1.861558

Проведем диагностику модели

# Assumption - mean of the residuals is = 0 
mean(model2$residuals)

## [1] 4.05783e-18

plot(model2$residuals)
abline(h=0,col='red')

Остатки распределены более меннее нормально.

Найдем выбросы.

outlierTest(model2)

##     rstudent unadjusted p-value Bonferroni p
## 413 4.270130         2.3419e-05     0.011850
## 372 4.052312         5.8859e-05     0.029782

У нас есть два выброса, которые можно удалить.

Boston_1 = Boston[-c(413,372),,drop=T]# drop the points
row.names(Boston_1)=1:nrow(Boston_1)

model3 = lm(log(medv) ~ crim + zn + chas + nox + rm + dis +
              ptratio + black + lstat,data=Boston_1)
par(mfrow=c(2,2))
plot(model3)

summary(model3)

## 
## Call:
## lm(formula = log(medv) ~ crim + zn + chas + nox + rm + dis + 
##     ptratio + black + lstat, data = Boston_1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.70878 -0.09807 -0.00993  0.08974  0.77161 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.6741127  0.1896496  19.373  < 2e-16 ***
## crim        -0.0084826  0.0011748  -7.220 1.97e-12 ***
## zn           0.0007490  0.0005193   1.442 0.149874    
## chas1        0.1199066  0.0335829   3.570 0.000391 ***
## nox         -0.6147940  0.1258689  -4.884 1.40e-06 ***
## rm           0.1123588  0.0157923   7.115 3.96e-12 ***
## dis         -0.0447458  0.0073176  -6.115 1.97e-09 ***
## ptratio     -0.0333834  0.0045429  -7.348 8.36e-13 ***
## black        0.0004409  0.0001050   4.197 3.20e-05 ***
## lstat       -0.0289046  0.0018797 -15.377  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1869 on 494 degrees of freedom
## Multiple R-squared:  0.7937, Adjusted R-squared:  0.7899 
## F-statistic: 211.1 on 9 and 494 DF,  p-value: < 2.2e-16

Новая модель:

R-squared составляет 0.7899
F-statistic составляет 211.1

После удаления выбросов, мы видим, что значение R-squared увеличилось. Однако один из параметров перестал быть значимым.

model3_step=stepAIC(model3,Boston_1,direction = 'both')

## Start:  AIC=-1680.97
## log(medv) ~ crim + zn + chas + nox + rm + dis + ptratio + black + 
##     lstat
## 
##           Df Sum of Sq    RSS     AIC
## <none>                 17.247 -1681.0
## - zn       1    0.0726 17.320 -1680.8
## - chas     1    0.4451 17.692 -1670.1
## - black    1    0.6151 17.862 -1665.3
## - nox      1    0.8329 18.080 -1659.2
## - dis      1    1.3054 18.552 -1646.2
## - rm       1    1.7673 19.014 -1633.8
## - crim     1    1.8202 19.067 -1632.4
## - ptratio  1    1.8853 19.132 -1630.7
## - lstat    1    8.2551 25.502 -1485.8

model3_step$anov

## Stepwise Model Path 
## Analysis of Deviance Table
## 
## Initial Model:
## log(medv) ~ crim + zn + chas + nox + rm + dis + ptratio + black + 
##     lstat
## 
## Final Model:
## log(medv) ~ crim + zn + chas + nox + rm + dis + ptratio + black + 
##     lstat
## 
## 
##   Step Df Deviance Resid. Df Resid. Dev       AIC
## 1                        494   17.24699 -1680.969

summary(model3_step)

## 
## Call:
## lm(formula = log(medv) ~ crim + zn + chas + nox + rm + dis + 
##     ptratio + black + lstat, data = Boston_1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.70878 -0.09807 -0.00993  0.08974  0.77161 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.6741127  0.1896496  19.373  < 2e-16 ***
## crim        -0.0084826  0.0011748  -7.220 1.97e-12 ***
## zn           0.0007490  0.0005193   1.442 0.149874    
## chas1        0.1199066  0.0335829   3.570 0.000391 ***
## nox         -0.6147940  0.1258689  -4.884 1.40e-06 ***
## rm           0.1123588  0.0157923   7.115 3.96e-12 ***
## dis         -0.0447458  0.0073176  -6.115 1.97e-09 ***
## ptratio     -0.0333834  0.0045429  -7.348 8.36e-13 ***
## black        0.0004409  0.0001050   4.197 3.20e-05 ***
## lstat       -0.0289046  0.0018797 -15.377  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1869 on 494 degrees of freedom
## Multiple R-squared:  0.7937, Adjusted R-squared:  0.7899 
## F-statistic: 211.1 on 9 and 494 DF,  p-value: < 2.2e-16

Новая модель:

R-squared составляет 0.7899
F-statistic составляет 211.1

par(mfrow=c(2,2))
plot(model3_step)

Мы видим, что нет разницы между model3 и model3_step, а также Р-квадрат имеет одинаковые значения. Теперь можем проанализировать выбросы модели.

outlierTest(model3)

##     rstudent unadjusted p-value Bonferroni p
## 369 4.303254         2.0307e-05     0.010235
## 372 4.266323         2.3833e-05     0.012012

Мы видим 2 выбросы, которые мы можем выбросить и протестировать модель еще раз.

Boston_2=Boston_1[-c(371,400,373),,drop=T]
row.names(Boston_2)=1:nrow(Boston_2)
model4 = lm(log(medv) ~ crim + zn + chas + nox + rm + dis +
                 tax + ptratio + black + lstat,data=Boston_2)
summary(model4)

## 
## Call:
## lm(formula = log(medv) ~ crim + zn + chas + nox + rm + dis + 
##     tax + ptratio + black + lstat, data = Boston_2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.70835 -0.10150 -0.00960  0.08772  0.79154 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.624e+00  1.899e-01  19.086  < 2e-16 ***
## crim        -7.880e-03  1.217e-03  -6.473 2.34e-10 ***
## zn           7.611e-04  5.218e-04   1.459  0.14530    
## chas1        1.062e-01  3.335e-02   3.185  0.00154 ** 
## nox         -5.614e-01  1.376e-01  -4.079 5.27e-05 ***
## rm           1.154e-01  1.550e-02   7.445 4.38e-13 ***
## dis         -4.340e-02  7.194e-03  -6.033 3.18e-09 ***
## tax         -6.026e-05  8.344e-05  -0.722  0.47055    
## ptratio     -3.255e-02  4.937e-03  -6.594 1.11e-10 ***
## black        4.396e-04  1.048e-04   4.195 3.23e-05 ***
## lstat       -2.858e-02  1.870e-03 -15.288  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1831 on 490 degrees of freedom
## Multiple R-squared:  0.7971, Adjusted R-squared:  0.7929 
## F-statistic: 192.4 on 10 and 490 DF,  p-value: < 2.2e-16

Новая модель:

R-squared составляет 0.7929
F-statistic составляет 192.4

par(mfrow=c(2,2))
plot(model4)

Мы видим, что мы снова увеличили мощность нашей модели, но потеряли значение по другой переменной. Мы можем запустить stepAIC и посмотреть, сможем ли мы его удалить.

outlierTest(model4)

##      rstudent unadjusted p-value Bonferroni p
## 369  4.527033         7.5217e-06    0.0037683
## 371  4.517479         7.8548e-06    0.0039352
## 372  3.982469         7.8571e-05    0.0393640
## 398 -3.953793         8.8281e-05    0.0442290

Мы видим 4 выбросы, которые мы можем выбросить и протестировать модель еще раз.

Boston_3 =Boston_2[-c(369,371,372,398),,drop=T]
row.names(Boston_3)=1:nrow(Boston_3)
Boston_3$chas <- as.factor(Boston_3$chas)
model5 = lm(log(medv) ~ crim + zn + chas + nox + rm + dis +
                 tax + ptratio + black + lstat,data=Boston_3)
summary(model5)

## 
## Call:
## lm(formula = log(medv) ~ crim + zn + chas + nox + rm + dis + 
##     tax + ptratio + black + lstat, data = Boston_3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.64945 -0.09613 -0.01079  0.08989  0.69041 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.421e+00  1.781e-01  19.207  < 2e-16 ***
## crim        -7.826e-03  1.132e-03  -6.915 1.48e-11 ***
## zn           5.723e-04  4.854e-04   1.179 0.239025    
## chas1        7.966e-02  3.144e-02   2.533 0.011610 *  
## nox         -5.005e-01  1.282e-01  -3.905 0.000108 ***
## rm           1.385e-01  1.468e-02   9.436  < 2e-16 ***
## dis         -3.810e-02  6.719e-03  -5.670 2.45e-08 ***
## tax         -1.143e-04  7.801e-05  -1.465 0.143473    
## ptratio     -3.225e-02  4.591e-03  -7.023 7.35e-12 ***
## black        4.324e-04  9.785e-05   4.418 1.23e-05 ***
## lstat       -2.677e-02  1.789e-03 -14.961  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1701 on 486 degrees of freedom
## Multiple R-squared:   0.82,  Adjusted R-squared:  0.8163 
## F-statistic: 221.4 on 10 and 486 DF,  p-value: < 2.2e-16

Новая модель:

R-squared составляет 0.8163
F-statistic составляет 221.4

par(mfrow=c(2,2))
plot(model5)

Мы улучшили еще немного нашу модель и теперь она объясняет 82% изменчивости. Нет предела совершенству, меня устраивает модель, объясняющая 82% изменчивости. Наибольший вклад вносят параметры: количество комнат, количество оксида азота, отношение учеников к учителям. Так что я бы посоветовала клиенту, строить большие дома рядом со школой и с парком неподалеку.

Финальный анализ модели.

residuals <- data.frame('Residuals' = model5$residuals)
res_hist <- ggplot(residuals, aes(x=Residuals)) + geom_histogram(color='black', fill='skyblue') + ggtitle('Histogram of Residuals')
res_hist

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Остатки распределены нормально.

plot(model5, col='Sky Blue')

# пересчитаем тестовый датасет
set.seed(123)
split <- sample.split(Boston,SplitRatio = 0.75) #assigns booleans to a new coloumn based on the split ratio
train <- subset(Boston,split==TRUE)
test <- subset(Boston,split==FALSE)

test$predicted.medv <- predict(model5,test)
pl1 <-test %>% 
  ggplot(aes(medv,predicted.medv)) +
  geom_point(alpha=0.5) + 
  stat_smooth(aes(colour='black')) +
  xlab('Значения medv') +
  ylab('Предсказанные значения medv') +
  theme_bw()

ggplotly(pl1)

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Проект 2

Подготовка к работе

2. Подготовим данные

3. Разведочный анализ.

Посчитаем корреляцию между признаками.

Визуализируем данные.

4. Построение модели, прогнозирующей стоимость дома

Задание 1. Построение полной модели.

Задание 2. Проведем диагностику полной модели.

Задание 3. Постройте график предсказаний стоимости от переменной, которая обладает

Дополнительное задание.

Финальный анализ модели.